{ "cells": [ { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import pandas\n", "import numpy\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will fit an ARX model of this form to a step response.\n", "\n", "$$y(k) = a_1 y(k - 1) + a_2 y(k - 2) + b_1 u(k - 1) + b_2 u(k - 2)$$" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ykuk
k
-10.0000
00.0001
10.0581
20.2171
30.3601
40.4881
50.6001
60.6921
70.7721
80.8331
90.8881
100.9251
\n", "
" ], "text/plain": [ " yk uk\n", "k \n", "-1 0.000 0\n", " 0 0.000 1\n", " 1 0.058 1\n", " 2 0.217 1\n", " 3 0.360 1\n", " 4 0.488 1\n", " 5 0.600 1\n", " 6 0.692 1\n", " 7 0.772 1\n", " 8 0.833 1\n", " 9 0.888 1\n", " 10 0.925 1" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = pandas.read_csv('../../assets/data.csv', index_col='k')\n", "data['uk'] = 1\n", "data.loc[0] = [0, 1] # input changes at t=0\n", "data.loc[-1] = [0, 0] # everything was steady at t < 0\n", "data = data.sort_index()\n", "data" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "y = data.yk\n", "u = data.uk" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXyV5Zn/8c9FCIR9SUiAJJiwE1G2AIpjRRHFuqBVXEvVYlGrxVFnOk5t3drfjK22tp2qLSpKFaGiVZDiVkdHqwIJi8giSQxb2EJOQlhDlnP//shipIEEOMlzznO+79fLFzk5T865osnXh/u5nus25xwiIuJfrbwuQEREmpeCXkTE5xT0IiI+p6AXEfE5Bb2IiM+19uqNExISXFpamldvLyISkZYvX17knOtxPF/jWdCnpaWRnZ3t1duLiEQkM9t8vF+jpRsREZ9T0IuI+JyCXkTE5xT0IiI+p6AXEfG5RoPezGaZWaGZrTnK82ZmvzezPDNbbWYjQ1+miIicqKac0b8ATDrG8xcBA2r+mQ48ffJliYhIqDTaR++c+8jM0o5xyGTgz6563vESM+tqZr2ccztCVGPoVJbDkqeg/IDXlYiItJhQ3DCVDGyt97ig5nP/FPRmNp3qs3769OkTgrc+Tvkfwt8frK2m5d9fRMQDoQj6hhKzwd1MnHMzgZkAmZmZLb/jSVFO9Z//ng8d4lv87UVETtrDx3+SGoqumwIgtd7jFGB7CF439AK50K67Ql5Eokoogn4h8L2a7pszgNKwXJ8HKMqFhAFeVyEi0qIaXboxs7nAeCDBzAqAB4FYAOfcH4HFwLeBPOAgcHNzFXvSinJhwAVeVyEi0qKa0nVzXSPPO+COkFXUXA7tgQOFOqMXkagTPXfGBvKq/1TQi0iUiZ6gL8qt/jNhoLd1iIi0sCgK+hxo1Rq6pXldiYhIi4qeoA/kQrd0iIn1uhIRkRYVPUGv1koRiVLREfTBKijOh/j+XlciItLioiPo92yGqnJdiBWRqBQdQV/XcaOlGxGJPqEYahb+1FopIhHKOcfW4kMs3Rgga1PxCb1GlAR9TvUws/bdva5EROSYgkFHbuF+lm0MsGxTCcs2Bti19zAAXdufWNdgdAR9IE/LNiISliqqgqzdvrc62DeWkL25mD0HKwDo2TmOsenxjEnvzpj07vTv0ZGYBxt5wQZER9AX5cJADTMTEe+VVVSxcsselm0sJmtTMcs3l3CoogqA9IQOXJjRsy7YU7q1w+zkN0nyf9DXDjOL1xm9iLS8vWUVLN9UwtKaYF9dsIeKKocZDO7ZmWtGpzI6rTuj07uR2CmuWWrwf9DXDTPThVgRaX679x0ma1MxyzZW/7N+516cg9gY47TkLkz7l76MSe/GqFO606Vdy9yp7/+gr90+UGv0IhJizjkKSg7VhXrWpmLyiw4A0C42hpGndOVfJwxkdHo3RqR2o12bGE/qjIKgz9UwMxEJmcJ9Zby7dlddsO8oLQOgS7tYRqd149ox1UsxQ5O7EBsTHrcqRUHQ52iYmYiclP2HK3lnzU7eWLWNT/KKCDpI7NSWMendGZvendHp3RmY2IlWrU7+wmlz8H/Qq7VSRE5ARVWQj3J288aq7by3bidlFUFSu7fjjnP7c9mw3vRP7BiSjpiW4O+gr6qsHmY28EKvKxGRCOCcY8WWPbyxcht/+2IHxQfK6dY+limjUrl8RG9G9ukWMeFen7+DvnaYmVorReQYvtq9nwUrt/HGqu1sKT5I29atmJiRxBUjkjl7QA/atA6PtfYT5e+gV2uliBxF4b4y3vx8BwtWbWN1QSmtDM7qn8BdEwZw4dCedGzrn3j0z3fSELVWikg9DV1UPS25Cz+9eAiXDetNYufmuWHJaz4P+lwNMxOJchVVQT7O3c3rK7++qJrSrfqi6uThyfRP7Oh1ic3O30EfyNOyjUgUqr2oumDVNhat/vqi6lWjUrhiRHLEXlQ9Uf4O+qIcddyIRJGjXVS9fHgy3xoY+RdVT5R/g/5QCRzYrTN6EZ9r6KLquH4JzJgwgAtPTaJTnG6W9G/QF9V03Ki1UsR3yiqqeGftTl5bsY1/5O4m6GBocmffX1Q9Uf4N+oC2DxTxm5xd+5i7bAuvr9zGnoMVpHRrxw/H9+fyEb3pn9jJ6/LCln+DviinZpjZKV5XIiIn4VB5FYtWb2de1laWby4hNsa48NSeXD+mD2f0jQ/b+TLhxMdBn6thZiIRbO32UuYt28obq7axr6ySvj06cP+3h/CdkcnEd2zrdXkRxb9Br9ZKkYiz/3Alb36+nXnLtvB5QSltWrfi4tN6cd2YPoxOi66WyFBqUtCb2STgd0AM8Kxz7tEjnu8DzAa61hxzn3NucYhrbbqqSgh8pdZKkQjgnOOLbaXMXbaFhau2c6C8ikFJnXjw0gyuGJFM1/ZtvC4x4jUa9GYWAzwJTAQKgCwzW+icW1fvsJ8CrzjnnjazDGAxkNYM9TbNns0QrNAZvUgY21tWwYKV25i7bCvrduylXWwMl5zei+vG9mFEaledvYdQU87oxwB5zrl8ADObB0wG6ge9AzrXfNwF2B7KIo9bUU3HjVorRcJK7R2rc5dtYdHq7ZRVBMno1ZmfXz6UycN701k9782iKUGfDGyt97gAGHvEMQ8B75rZj4AOwPkNvZCZTQemA/Tp0+d4a226utZKBb1IONhzsJy/rtjGvKwt5OzaT4c2MVwxIoXrxqRyWnIXnb03s6YEfUP/BdwRj68DXnDO/drMzgReNLOhzrngN77IuZnATIDMzMwjXyN0inKhfbyGmYl4yDnHso3FzF22hcVrdlJeGWRYalce/c5pXDqsNx18NAY43DXl33QBkFrvcQr/vDQzDZgE4Jz7zMzigASgMBRFHreiXC3biHgksP8wf12xjblZW8jffYBOca25dnQq147uQ0bvzo2/gIRcU4I+CxhgZunANuBa4PojjtkCTABeMLMhQBywO5SFHpdALgyc5Nnbi0SbYNDxWX6Al5dt4d21O6mocmSe0o0fTunPxaf1ol2bGK9LjGqNBr1zrtLM7gTeobp1cpZzbq2ZPQJkO+cWAvcCz5jZ3VQv69zknGu+pZljqRtmpjN6kea2t6yCOUu2MHfZFrYUH6RLu1imnpHGtWNSGZikkQThokmLZDU98YuP+NwD9T5eB5wV2tJOUJG2DxRpbofKq5j92Sb++H9fsedgBWPTu3PvBQO58NSexMXq7D3c+O9qSO32gVqjFwm5w5VVzFu2lT98kMfufYcZP6gH904cxGkpXbwuTY7Bf0EfyNUwM5EQq6wK8tqKAn7/fh7b9hxibHp3nrphJKPT1NkWCfwX9EW50L2vhpmJhEAw6Fj0xQ6eeC+HjUUHGJbShUevPI1/6Z+g3vcI4s+g17KNyElxzvH39YX8+t0NfLlzH4OSOjFz6igmZiQp4COQv4K+qhKK82HQRV5XIhKRnHN8khfgsXc38PnWPaQndOB31w7n0tN7a+57BPNX0NcNM9MZvcjxWr65mMfe2cCS/GJ6d4njl1eexpUjU2gdE50bavuJv4K+SNsHihyvNdtK+fW7G/hgw24SOrbloUszuG5sH9q2VpukX/gr6GuHmcX397YOkQiQV7iP37yXw+IvdtKlXSz/MWkwN447hfZt/BUL4regL8rRMDORRmwJHOS37+fwxspttIuNYcaEAdxydrpGBPuYz4Je2weKHM3O0jL+539z+UvWVmJaGbec3ZfbzulH9w7awcnvfBb0Oeq4ETlCYP9hnv7wK15cspmgc1w3pg93ntefpM5xXpcmLcQ/QX+wGA4W6YxepEbpoQqe/TifWf/YyKGKKr4zMoW7JgwgtXt7r0uTFuafoA/UDjNTa6VEt4PllTz/ySZmfpRP6aEKLj69F3efP5D+iR29Lk084p+g1z6xEuXKKqp4eekWnvowj6L95UwYnMg9Fwzk1N4aOBbt/BP0gVxoFathZhJ1gkHH6yu38fi7G9hRWsa4fvH8aeogRp3SzevSJEz4J+iLcqF7uoaZSVRZs62UBxeuZfnmEoalduXXU4Yxrn+C12VJmPFX0OtCrESJ0oMVPP7uBuYs3Uy39m147KrTuXJkiubRSIP8EfQaZiZRIhh0zF++lV++vYE9B8v53plp3D1xIF3a6W+ycnT+CPq6YWY6oxf/Wl2whwcWrGXV1j2MTuvGw5eNJaN3Z6/Lkgjgj6CvG2amjhvxn5ID5Tz27gbmLttCfIe2/ObqYVwxIllz4aXJfBL0tfvEapiZ+EdV0PGXrK386p0v2VdWyc3j0vnXiQM0k0aOmz+CPpAL7RM0zEx8Y9XWPTywYA2rC0oZk96dn08eyqCenbwuSyKUP4K+KFfLNuILgf2HeeydDfwleys9Orbld9cO57JhvbVMIyfFP0E/+NteVyFywqqCjpeXbubxd3M4cLiSH5zdlxkTBtCxrT9+RcVbkf9TVDvMTKMPJEIt31zCAwvWsHb7Xsb1i+fhy05lQJKWaSR0Ij/oNcxMIlTR/sM8+taXvLq8gJ6d4/jD9SO4+LReWqaRkIv8oNc+sRJhKquCvLRkM79+L4eyiipuO6cfPzqvPx20TCPNJPJ/sopyqoeZddUwMwl/WZuK+dkba/hy5z7OHpDAQ5edSr8eGh8szSvygz6QB937QkzkfyviX4X7ynh08Zf8deU2eneJ4+kbRjJpaE8t00iLiPx0LMrRso2ErcqqILM/28xv38vhcGWQO87txx3n9qd9m8j/1ZPI0aSfNjObBPwOiAGedc492sAxVwMPAQ743Dl3fQjrbFhVJRRvhMEXN/tbiRyvJfkBHlywlg279nHOwB48dNmppCd08LosiUKNBr2ZxQBPAhOBAiDLzBY659bVO2YA8J/AWc65EjNLbK6Cv6F2mJlaKyWM7Npbxn8tXs+CVdtJ7tqOmVNHMTEjScs04pmmnNGPAfKcc/kAZjYPmAysq3fMD4AnnXMlAM65wlAX2qDaGTdqrZQwUBV0PP/JRp54L4eKoGPGhAHcfk4/2rWJ8bo0iXJNCfpkYGu9xwXA2COOGQhgZp9QvbzzkHPu7SNfyMymA9MB+vTpcyL1flPdPrEaZibe2llaxt1/WcVn+QHOG5zIg5dmcEq8lmkkPDQl6Bv6+6Zr4HUGAOOBFOBjMxvqnNvzjS9ybiYwEyAzM/PI1zh+RTkaZiaee3/9Lv5t/ueUVQR57KrTuWpUipZpJKw0JegLgNR6j1OA7Q0cs8Q5VwFsNLMNVAd/VkiqPJpAnjpuxDNlFVU8+taXvPDpJjJ6deZ/rh+hnngJS62acEwWMMDM0s2sDXAtsPCIY94AzgUwswSql3LyQ1log4pyIEHLNtLy8gr3c8VTn/LCp5v4/lnpvH7HOIW8hK1Gz+idc5VmdifwDtXr77Occ2vN7BEg2zm3sOa5C8xsHVAF/LtzLtCchVcPMwuo40ZalHOOV7K38tDCdbRrE8OsmzI5b3CS12WJHFOT+uidc4uBxUd87oF6Hzvgnpp/WkbdMDMt3UjL2FtWwU/++gWLVu/grP7xPHH1cBI7x3ldlkijIvf2PLVWSgtasaWEGXNXsqO0jB9PGsRt3+pHq1a64CqRIYKDPlfDzKTZVQUdf/y/r/jNezn07hrH/NvOZGSfbl6XJXJcIjvoNcxMmtGuvdW98Z9+FeDSYb35f1cM1cbcEpEiNyUDuVqfl2aj3njxk8gM+qoKDTOTZnG4sor/XqzeePGXyAz6Eg0zk9DLK9zPjLkrWbdjLzeflcZ9Fw2mbWvNqZHIF5lBH9D2gRI6zjnmZxfw4MK1tGsTw3M3ZjJhiHrjxT8iM+jrWit1V6ycnPq98eP6xfPENcNJUm+8+EyEBn0udOgB7dTmJifuyN74W7/Vjxj1xosPRW7Qa31eTlD93vheXdQbL/4XmUEfyFXHjZyQ+r3xl5zei//6zmnqjRffi7ygrx1mpguxcpzq98b/6qrTmaLeeIkSkRf0dbtKaelGmka98RLtIi/o61orFfTSOPXGi0Ri0BflaJiZNEq98SJfi8Cgz4P4fhpmJkel3niRb4q8tNQwMzmGTUUHmDY7i02Bg/z7hYO47Rz1xotEVtBXVUBxvlorpUGf5hVx+5wVtDKYc8tYzugb73VJImEhsoK+ZDMEK3VGL//kpSWbeWjhWtITOvDcjaPpE9/e65JEwkZkBX3tjBu1VkqNyqogP1+0jtmfbebcQT34/XUj6KQboES+IbKCvq61UsPMBEoPVnDHyyv4R14RPzg7nfsuGqL1eJEGRFbQF+VomJkAkL97P7fMzmZryUF+deXpXD061euSRMJWhAV9npZthH/kFvHDOctpHdOKObecwZj07l6XJBLWWnldwHEJ5OqO2Cj34mebuPH5ZfTq0o4Fd5ylkBdpgsg5o68bZqagj0YVVUEefnMtLy3ZwvlDEvnttSPo2DZyfnxFvBQ5vylF2j4wWu05WM4dL6/gk7wAt57Tlx9fOFgXXUWOQwQFfW1rpTpuokle4X5umZ3F9j1lPD5lGFeNSvG6JJGIEzlBH8iFmDYaZhZFPsrZzR0vr6Bt61bMnT6WUadoPV7kRERO0BflQfe+GmYWBZxzvPDpJn6+aB0Dkzrx7I2ZpHTTna4iJypyUrMoB3oM8roKaWYVVUEeWLCWucu2cEFGEk9cM5wOuugqclIi4zeoqgJKNsKQS72uRJpRyYFybp+znCX5xfxwfD/+7YJBtNJFV5GT1qQ+ejObZGYbzCzPzO47xnFXmZkzs8zQlQiUbKoZZqbWSr/K3bWPyU9+woote/jtNcP58aTBCnmREGn0jN7MYoAngYlAAZBlZgudc+uOOK4TMANYGvIq1Vrpax9sKGTGyytpGxvDvOlnMLKPRlyIhFJTzujHAHnOuXznXDkwD5jcwHE/B34FlIWwvmpqrfQl5xzPfpzPtBeySO3engV3nqWQF2kGTQn6ZGBrvccFNZ+rY2YjgFTn3KJjvZCZTTezbDPL3r17d9OrDOTWDDPr2vSvkbBWXhnkvte+4Bd/W8/EjCRevf1Mkru287osEV9qysXYhhZKXd2TZq2AJ4CbGnsh59xMYCZAZmama+TwrxXladnGRwL7D3P7SytYtqmYH53Xn7vPH6j1eJFm1JSgLwDqz4BNAbbXe9wJGAp8aGYAPYGFZnaZcy47JFUW5ajjxic27NzHtNlZFO47zO+uHc7k4cmNf5GInJSmBH0WMMDM0oFtwLXA9bVPOudKgYTax2b2IfBvIQv5AwE4VKwzeh94f/0uZsxdSfu2rXnl1jMZnqqlOJGW0GjQO+cqzexO4B0gBpjlnFtrZo8A2c65hc1aYd2uUmqtjFTOOZ75OJ//futLTu3dmWe+l0mvLlqPF2kpTbphyjm3GFh8xOceOMqx40++rHqKFPSR7HBlFfe/voZXlxdw8Wm9eHzKMNq1ifG6LJGoEv53xmqYWcQq2n+Y215cTvbmEu6aMIC7JgzQRVcRD4R/0BflVg8za6WzwEiyqegA331uKbv3HeYP14/gktN7e12SSNSKjKBPHOx1FXIc1m3fy/dmLSPoHK/ceibDdNFVxFPhvWds7TAzbQgeMbI3FXPNzM+IjTGFvEiYCO8z+rphZmqtjAQfbCjk9peW07tLO168ZazudBUJE+Ed9LUzbtRxE/YWfr6de/6yikE9OzH7+2NI6NjW65JEpEaYB31Na6WGmYW1l5Zs5mcL1jA6rTvP3phJ57hYr0sSkXrCO+gDudAhUcPMwpRzjqc+/IrH3tnAhMGJPHnDSOJi1R0lEm7CO+iLcrVsE6acc/zX4vU88/FGLh/em8emDCM2Jryv7YtEq/D+zVTQh6XKqiD/8dpqnvl4IzeeeQq/uXq4Ql4kjIXvGX3tMDO1VoaVw5VV3DV3FW+v3cmMCQO4+/wB1EwtFZEwFb5BH9D2geFm/+FKbn0xm0/yAjxwSQbf/5d0r0sSkSYI36CvG2amjptwUHKgnJteyGLNtlJ+PWUYV45K8bokEWmiMA76HA0zCxM7S8uY+txSNhcf5I/fHcXEjCSvSxKR4xC+QR/Ig+79NMzMY7XDyUoOlPPCzaMZ1y+h8S8SkbASvkFflAOJQ7yuIqrVDierCgaZO/0MTk/R/QwikSg8e+KqKqrn3OhCrGfqDyebf9uZCnmRCBaeZ/S1w8zUWumJDzcUcttLy+nVpR0vThtDSrf2XpckIichPINew8w88+bn27n7L6sYmNSJP0/TcDIRPwjToNcwMy/MWbqZn76xhtGndOfZmzScTMQvwjfoNcysxdQfTnbe4ESevH6kNvAW8ZHwDPpAri7EthDnHP/91pfM/CifycN787iGk4n4Tnj+Rhfl6I7YFlA7nGzmR/l878xTeELDyUR8KfzO6A8E4FCJzuib2TeGk53Xn7snDtRwMhGfCr+grx1mptbKZnPgcCXTa4aT/eySDKZpOJmIr4Vf0Ne1VmrppjmUHCjn5hey+GJbKY9PGcZVGk4m4nthGPS5GmbWTOoPJ3v6hpFccGpPr0sSkRYQnkGvYWYhp+FkItEr/II+kAuJGV5X4Svrd+xl6nPVw8le/sEZDEvV/Qki0SS8eunqhpnpQmyorNq6h2v+9BmtW1UPJ1PIi0SfJgW9mU0ysw1mlmdm9zXw/D1mts7MVpvZ+2Z2YgvsxRs1zCyE1mwrZepzS+navg2v3n4m/RM7eV2SiHig0aA3sxjgSeAiIAO4zsyOXFtZCWQ6504HXgV+dULVaJ/YkPly516mPreUznGxvPyDsZpAKRLFmnJGPwbIc87lO+fKgXnA5PoHOOc+cM4drHm4BDixnj21VoZEXuF+vvvsUtq0bqWQF5EmBX0ysLXe44Kazx3NNOCthp4ws+lmlm1m2bt37/7nA4ryoGMSxHVpQlnSkE1FB7j+mSWA8fIPzuCU+A5elyQiHmtK0Dd0X7xr8ECz7wKZwGMNPe+cm+mcy3TOZfbo0eOfDyjK0fr8SSgoOcgNzy6loirInFvG0q9HR69LEpEw0JSgLwBS6z1OAbYfeZCZnQ/cD1zmnDt8QtUEctVxc4J2lB7i+meWsq+sghenjWVQT114FZFqTQn6LGCAmaWbWRvgWmBh/QPMbATwJ6pDvvCEKqkbZqagP16F+8q44ZmlFB8o58/TxjI0WUtfIvK1RoPeOVcJ3Am8A6wHXnHOrTWzR8zssprDHgM6AvPNbJWZLTzKyx1d7YVYLd0cl+ID5Xz32aXsKC3j+ZtHM1x98iJyhCbdGeucWwwsPuJzD9T7+PyTrqSutVJB31SlByv47rNL2Rw4yPM3jWZ0WnevSxKRMBQ+d8YW5UBMW+jax+tKIsK+sgq+N2speYX7+dPUUYzrr9k1ItKwMAr6PIjXMLOmOHC4kpufz2Lt9r08ecNIxg9K9LokEQlj4RP0gVyI141SjTlUXsW02Vms2FLC768bwcSMJK9LEpEwFx5BX1lePedG6/PHVFZRxfQXs1m6sZjfXD2cb5/Wy+uSRCQChEfQl2wCV6UZN8dQXhnkzpdX8HFuEb/8zulcPuJYNyeLiHwtPIJerZXHVFkV5K55K/n7+kJ+fvlQrh6d2vgXiYjUCI+NR+paK7VGf6SqoOPe+Z/z1pqd/PTiIUw9Q1ssijRVRUUFBQUFlJWVeV3KcYuLiyMlJYXY2NiTfq3wCPqiXA0za0Aw6LjvtdUsWLWdH08axC1n9/W6JJGIUlBQQKdOnUhLS8OsobFd4ck5RyAQoKCggPT09JN+vTBZusnV+vwRnHM8sHAN85cXMGPCAH44Xn/bETleZWVlxMfHR1TIA5gZ8fHxIfubiPdB71zN1EoFWS3nHD9ftJ6Xlmzh1nP6cvf5unYhcqIiLeRrhbJu74P+YADK9qi1soZzjl+9s4FZn2zkpnFp3DdpcMT+oIpIePA+6Iu0fWB9v38/j6c//Irrx/bhwUszFPIiPvPhhx9yySWXtOh7hkHQ17ZWaunm6Q+/4om/53DVqBR+MXmoQl5EQsL7rptAroaZAbP+sZFfvv0llw7rzS+vPJ1WrRTyIqH08JtrWbd9b0hfM6N3Zx689NSjPv+zn/2MhIQE7rrrLgDuv/9+kpK+HluSlZXF9OnTee211+jbt/m66sLgjF7DzF5asplHFq1j0qk9+c3Vw4hRyIv4wrRp05g9ezYAwWCQefPmkZxcfVf7p59+ym233caCBQuaNeQhHM7oi3Kg51Cvq/DMK9lb+ekbazhvcCK/v24EsTHe/79XxI+OdebdXNLS0oiPj2flypXs2rWLESNGEB8fz/r165k+fTrvvvsuvXv3bvY6vA36yvLqOTenXuFpGV5ZsGob//Haas4ekMBTN4ykTWuFvIjf3HLLLbzwwgvs3LmT73//+wD06tWLsrIyVq5c2SJB722ylGysGWYWfa2Vb6/ZwT2vfM6YtO7MnJpJXGz0Ll2J+NkVV1zB22+/TVZWFhdeeCEAXbt25W9/+xs/+clP+PDDD5u9Bm/P6Iuic/vA99fv4kdzVzI8tSuzbhpNuzYKeRG/atOmDeeeey5du3YlJubr3/WkpCTefPNNLrroImbNmsXYsWObrQaPgz76plZ+lLOb219awZBenXn+5tF0aOv9ZRIRaT7BYJAlS5Ywf/58AMaPH8/48eMB6NOnD2vXrm32GrxdugnkQceeENfZ0zJaymdfBZj+Yjb9Ejvy5++PoXPcyU+lE5HwtW7dOvr378+ECRMYMMC7E1rvl26iZNnm07wibvlzNqnd2vPStDF0bd/G65JEpJllZGSQn5/vdRleB32Orztuig+Us2DVNuZnF7Bux17SEzow55axxHds63VpIhJFvAv6YGXNMDN/zbipqAryfxt28+ryAt7/chcVVY7TU7rwyORTuXxEspZrRKTFeRf0lTVzln2ydLNh5z5eXb6V11dup2j/YRI6tuGmcWlcOSqFwT2j4xqEiIQnD4P+cPWfERz0ew6Ws/Dz7by6vIDVBaW0bmVMGJLIlFGpnDOoh+5yFZGw4O0ZfUxb6BJZG11XVgX5OLeIV5cX8N66XZRXBcno1ZkHLslg8vDeWn8XkSYZP348jz/+OJmZmc3+Xt6e0ccPjiZjMRoAAAdoSURBVJhhZnmF+5i/vIDXV2yjcN9hundoww1n9OGqUSmc2lt73YpI+PL2jD4hvGfQlx6qYNHq7czPLmDV1j3EtDLOHZTIVaNSOG9wombTiESSt+6DnV+E9jV7ngYXPXrMQzZt2sQll1zCmjVrAHj88cfZv39/3fPBYJCbb76Z1NRUfvGLX4S2vhoeBn15WN4RWxV0fJJXvTTzztqdHK4MMiipEz+9eAiThyfTo5OWZkQkNCorK7nhhhsYOnQo999/f7O9j4d99C6sWivzd+/ntRUF/HXFNnaUltGlXSzXjE5lyqhUhiZ31m5PIpGukTNvL9x6661cffXVzRry0MSgN7NJwO+AGOBZ59yjRzzfFvgzMAoIANc45zY1+sIeL93sK6vgb6t38OryArI3l9DK4JyBPfjZJRlMGJJI29aRcf1ARMJX69atCQaDdY/LysrqPh43bhwffPAB9957L3Fxcc1XQ2MHmFkM8CQwESgAssxsoXNuXb3DpgElzrn+ZnYt8Evgmkbf3YOlm2DQsSQ/wPzlBby1ZgdlFUH6J3bkvosG850RySR2br5/2SISfZKSkigsLCQQCNCxY0cWLVrEpEmTgOodqD766COmTJnC66+/TuvWzbPI0pRXHQPkOefyAcxsHjAZqB/0k4GHaj5+FfiDmZlzzh3tRStpzcSnVp5Q0Sdjz6EKdu87TKe41lw5MoUpmakMS+mipRkRaRaxsbE88MADjB07lvT0dAYPHvyN5++55x5KS0uZOnUqc+bMoVWr0Dd5NCXok4Gt9R4XAEcOTq47xjlXaWalQDxQVP8gM5sOTAfo26sb45I6nmDZJ65NTCvOG5LEBRlJ2uxDRFrEjBkzmDFjxlGff/jhh5v1/ZsS9A2d6h55pt6UY3DOzQRmAmRmZrqnbhjVhLcXEZGT0ZS/IxQA9W9fTQG2H+0YM2sNdAGKQ1GgiIicnKYEfRYwwMzSzawNcC2w8IhjFgI31nx8FfC/x1qfFxFpKZEaRaGsu9Ggd85VAncC7wDrgVecc2vN7BEzu6zmsOeAeDPLA+4B7gtZhSIiJyguLo5AIBBxYe+cIxAIhKzl0rz6F5CZmemys7M9eW8RiQ4VFRUUFBR8o3c9UsTFxZGSkkJs7Df3sDCz5c6545qEpp2pRcS3YmNjSU9P97oMz2kql4iIzynoRUR8TkEvIuJznl2MNbPdwGZP3hwSOOKuXZ+Ltu8X9D1Hi2j8ngc55zodzxd4djHWOdfDq/c2s+zjvWodyaLt+wV9z9EiWr/n4/0aLd2IiPicgl5ExOeiNehnel1AC4u27xf0PUcLfc9N4NnFWBERaRnRekYvIhI1FPQiIj4XVUFvZpPMbIOZ5ZmZ7ydsmlmqmX1gZuvNbK2Z3eV1TS3FzGLMbKWZLfK6lpZgZl3N7FUz+7Lmv/eZXtfU3Mzs7pqf6zVmNtfMfLfhs5nNMrNCM1tT73Pdzew9M8ut+bNbY68TNUFfb5Pzi4AM4Dozy/C2qmZXCdzrnBsCnAHcEQXfc627qB6rHS1+B7ztnBsMDMPn37uZJQMzgEzn3FAghuq9MvzmBWDSEZ+7D3jfOTcAeJ8mjIWPmqCn3ibnzrlyoHaTc99yzu1wzq2o+Xgf1b/8yd5W1fzMLAW4GHjW61pagpl1Br5F9b4QOOfKnXN7vK2qRbQG2tXsateef975LuI55z7in3frmwzMrvl4NnB5Y68TTUHf0Cbnvg+9WmaWBowAlnpbSYv4LfBjIOh1IS2kL7AbeL5muepZM+vgdVHNyTm3DXgc2ALsAEqdc+96W1WLSXLO7YDqkzkgsbEviKagb9IG5n5kZh2B14B/dc7t9bqe5mRmlwCFzrnlXtfSgloDI4GnnXMjgAP4fJe3mnXpyUA60BvoYGbf9baq8BVNQd+UTc59x8xiqQ75Oc65v3pdTws4C7jMzDZRvTx3npm95G1Jza4AKHDO1f5t7VWqg9/Pzgc2Oud2O+cqgL8C4zyuqaXsMrNeADV/Fjb2BdEU9E3Z5NxXzMyoXrdd75z7jdf1tATn3H8651Kcc2lU/zf+X+ecr8/0nHM7ga1mNqjmUxOAdR6W1BK2AGeYWfuan/MJ+PwCdD0LgRtrPr4RWNDYF0TNVoLOuUozq93kPAaY5Zxb63FZze0sYCrwhZmtqvncT5xziz2sSZrHj4A5NScx+cDNHtfTrJxzS83sVWAF1d1lK/HhOAQzmwuMBxLMrAB4EHgUeMXMplH9P7wpjb6ORiCIiPhbNC3diIhEJQW9iIjPKehFRHxOQS8i4nMKehERn1PQi9Qws7T6UwJF/EJBLyLicwp6kQaYWd+aAWGjva5F5GQp6EWOUDNK4DXgZudcltf1iJysqBmBINJEPaieHXJlFIzIkCihM3qRbyqlet+Cs7wuRCRUdEYv8k3lVO/Y846Z7XfOvex1QSInS0EvcgTn3IGaDUzeM7MDzrlGx8CKhDNNrxQR8Tmt0YuI+JyCXkTE5xT0IiI+p6AXEfE5Bb2IiM8p6EVEfE5BLyLic/8feKtLKYGBYIsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We effectively have the following equations (I repeat the equation here for convenience):\n", "\n", "$$y(k) = a_1 y(k - 1) + a_2 y(k - 2) + b_1 u(k - 1) + b_2 u(k - 2)$$" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.058 = a1×0.0 + a2×0.0 + b1×1 + b2×0\n", "0.22 = a1×0.058 + a2×0.0 + b1×1 + b2×1\n", "0.36 = a1×0.22 + a2×0.058 + b1×1 + b2×1\n", "0.49 = a1×0.36 + a2×0.22 + b1×1 + b2×1\n", "0.6 = a1×0.49 + a2×0.36 + b1×1 + b2×1\n", "0.69 = a1×0.6 + a2×0.49 + b1×1 + b2×1\n", "0.77 = a1×0.69 + a2×0.6 + b1×1 + b2×1\n", "0.83 = a1×0.77 + a2×0.69 + b1×1 + b2×1\n", "0.89 = a1×0.83 + a2×0.77 + b1×1 + b2×1\n", "0.93 = a1×0.89 + a2×0.83 + b1×1 + b2×1\n" ] } ], "source": [ "for k in range(1, 11):\n", " print(f'{y[k]:.2} = a1×{y[k - 1]:.2} + a2×{y[k - 2]:.2} + b1×{u[k - 1]} + b2×{u[k - 2]}')" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ky[k]y[k-1]y[k-2]u[k-1]u[k-2]
010.0580.0000.00010
120.2170.0580.00011
230.3600.2170.05811
340.4880.3600.21711
450.6000.4880.36011
560.6920.6000.48811
670.7720.6920.60011
780.8330.7720.69211
890.8880.8330.77211
9100.9250.8880.83311
\n", "
" ], "text/plain": [ " k y[k] y[k-1] y[k-2] u[k-1] u[k-2]\n", "0 1 0.058 0.000 0.000 1 0\n", "1 2 0.217 0.058 0.000 1 1\n", "2 3 0.360 0.217 0.058 1 1\n", "3 4 0.488 0.360 0.217 1 1\n", "4 5 0.600 0.488 0.360 1 1\n", "5 6 0.692 0.600 0.488 1 1\n", "6 7 0.772 0.692 0.600 1 1\n", "7 8 0.833 0.772 0.692 1 1\n", "8 9 0.888 0.833 0.772 1 1\n", "9 10 0.925 0.888 0.833 1 1" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pandas.DataFrame([(k, \n", " y[k], \n", " y[k-1], \n", " y[k-2], \n", " u[k-1], \n", " u[k-2]) for k in range(1, 11)], columns=['k', 'y[k]', 'y[k-1]', 'y[k-2]', 'u[k-1]', 'u[k-2]'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We notice that these equations are linear in the coefficients. We can define $\\beta= [a_1, a_2, b_1, b_2]^T$. Now, to write the above equations in matrix form $Y = X \\beta $, we define" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.058],\n", " [0.217],\n", " [0.36 ],\n", " [0.488],\n", " [0.6 ],\n", " [0.692],\n", " [0.772],\n", " [0.833],\n", " [0.888],\n", " [0.925]])" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Y = numpy.atleast_2d(y.loc[1:]).T\n", "Y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To build the coefficient matrix we observe that there are two blocks of constant diagonals (the part with the $y$s and the part with the $u$s). Matrices with constant diagonals are called Toeplitz matrices and can be constructed with the `scipy.linalg.toeplitz` function." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "import scipy\n", "import scipy.linalg" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0. , 0. ],\n", " [0.058, 0. ],\n", " [0.217, 0.058],\n", " [0.36 , 0.217],\n", " [0.488, 0.36 ],\n", " [0.6 , 0.488],\n", " [0.692, 0.6 ],\n", " [0.772, 0.692],\n", " [0.833, 0.772],\n", " [0.888, 0.833]])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X1 = scipy.linalg.toeplitz(y.loc[0:9], [0, 0])\n", "X1" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 0],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1],\n", " [1, 1]])" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X2 = scipy.linalg.toeplitz(u.loc[0:9], [0, 0])\n", "X2" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "X = numpy.hstack([X1, X2])" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0. , 0. , 1. , 0. ],\n", " [0.058, 0. , 1. , 1. ],\n", " [0.217, 0.058, 1. , 1. ],\n", " [0.36 , 0.217, 1. , 1. ],\n", " [0.488, 0.36 , 1. , 1. ],\n", " [0.6 , 0.488, 1. , 1. ],\n", " [0.692, 0.6 , 1. , 1. ],\n", " [0.772, 0.692, 1. , 1. ],\n", " [0.833, 0.772, 1. , 1. ],\n", " [0.888, 0.833, 1. , 1. ]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another option is to use the loop from before to construct the matrices. This is a little more legible but slower:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "Y = []\n", "X = []\n", "for k in range(1, 11):\n", " Y.append([y[k]])\n", " X.append([y[k - 1], y[k - 2], u[k - 1], u[k - 2]])\n", "Y = numpy.array(Y)\n", "X = numpy.array(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We solve for $\\beta$ as we did for linear regression:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0.98464753],\n", " [-0.12211256],\n", " [ 0.058 ],\n", " [ 0.10124916]])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta, _, _, _ = numpy.linalg.lstsq(X, Y, rcond=None)\n", "beta" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }